Vertex dynamics in finite two dimensional square spin ices 
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Local magnetic ordering in artificial spin ices is discussed from the point of view of how geometrical 
frustration controls dynamics and the approach to steady state. We discuss the possibility of using a 
particle picture based on vertex configurations to interpret time evolution of magnetic configurations. 
Analysis of possible vertex processes allows us to anticipate different behaviors for open and closed 
edges and the existence of different field regimes. Numerical simulations confirm these results and 
also demonstrate the importance of correlations and long range interactions in understanding particle 
population evolution. We also show that a mean field model of vertex dynamics gives important 
insights into finite size effects. 

PACS numbers: 75.50.Lk, 75.75.-c, 75.78.-n 



o 

a 

i 



c 

o 

(N 
> 

in 

(N 

\o 

(N 
O 

o 



X 



Introduction. — Spin ices are geometrically frustrated 
magnetic systems, where interaction energies are mini- 
mized by local arrangements of spins resembling the ice 
rule for the ground state of solid water, i.e., two- in, two- 
out spin configurations [2, 0] . Spin ices display several 
interesting features including zero point entropy [s| , spin 
freezing, and hysteresis [3] • The role of long range inter- 
actions in three dimensional spin ices is also interesting, 
as the long range dipolar interactions between spins lead 
to a state that can be described using the short range ice 
rules @- 

There is now a growing interest in two dimensional ar- 
tificial spin ices, which consist of finite arrays of elongated 
magnetic dots whose magnetizations can be well approx- 
imated by Ising spins [a, @|. Artificial spin ices have 
several differences with their three dimensional counter- 
parts. For example, their geometry can be controlled ex- 
perimentally and their magnetization configurations 
can be imaged directly using scanning probe techniques. 
Furthermore, artificial spin ices allow the study of frus- 
tration at room temperature because of the thermal sta- 
bility of the relatively large dots. In fact, the anisotropy 
barrier preventing the macrospin associated with a dot 
from flipping is much larger than room temperature, so 
dynamics can only be induced by applying a magnetic 
field H = hO. We can expect that a small h does not af- 
fect the configuration, while a large h dominates dipolar 
interactions and spins simply follow the field Inter- 
esting collective effects arise when \h — h c \ ~ g, where h c 
is the minimum field required to flip an isolated spin and 
g is of the order of the dipolar interactions. 

Until now, studies of artificial spin ices have focused on 
the best procedures to attain low energy states [s| or the 
properties of 'final' states of demagnetization protocols 
{ioi \Ti\ . Very little is known about the time-dependent 
dynamical evolution of spin ices in response to a magnetic 
field, either experimentally or theoretically. In this Letter 
we show that the time evolution of square artificial spin 



ices can be described in terms of vertex configurations, 
and that these interact in a way that in principle can be 
predicted and understood. 

We do this by first introducing the essentials of a vertex 
population model, which predicts some general features 
of the dynamics. We show that a mean field approxi- 
mation allows quantitative evaluation of important array 
size effects. However, correlations, the long range part of 
dipolar interactions, and edge effects can be nonnegligi- 
ble and are studied carefully using numerical simulations. 
In all of our models, a rotating field of constant strength 
is applied to finite arrays and we discuss how this leads 
to effective demagnetization if h is correctly tuned. 

Vertex population model. — In this work, we consider 
only square artificial spin ices. Unlike the spins on 
the tetrahedra of three dimensional ices, the four spins 
around a vertex of a square ice are not equivalent (Fig. 
Hk). As in Ref. 0, we classify the vertex configurations 
into four types of increasing energy Ei (i — 1, . . . , 4) (Fig. 
HJd). Type I and 2 vertices obey the two- in, two-out ice 
rule but E\ < E%, because type 1 vertices minimize the 
energy of stronger interactions. 

The arrays can have open or closed edges (Fig. [Ib,d). 
Unlike in the open edge case, the configuration of an 
array with closed edges cannot be fully described by the 
configuration of its 4-island vertices and three types of 
3-island edge vertices must be introduced (Fig. [JJ). 

We now discuss the essentials of a vertex population 
model for open edge arrays. After sample saturation in 
the x direction, a field of constant modulus h is applied 
and rotated anticlockwise (Fig. QJ). The system evolves 
via single spin flips rather than so-called 'loop moves' 12 1 
so that ice-rule disobeying vertices can arise, a necessary 
condition to describe dynamics. Let us suppose that the 
condition for flipping a spin is —H\ ot ■ > h c . if* ot is 
the total field acting on Si, and is the sum of the external 
field H and the dipolar field Hf ip due to all other spins. 
We work in reduced units where the energy of a nearest- 
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FIG. 1: The geometry of square artificial spin ice. (a) The 
nearest and next nearest neighbor interactions at each 4- 
island vertex of the array; (b) examples of each type of 4- 
island vertex; (c) 'open' edge geometry which can be described 
fully by the 4-island vertex configuration; (d) 'closed' edge ge- 
ometry with vertices containing three islands; (e) examples of 
each type of 3-island vertex; (f) the magnitude and direction 
of the applied magnetic field; (g) the single spin flip one- and 
two-vertex processes that are energetically and topologically 
allowed. The spin that flips is black, while spins in type 1, 2, 
and 3 vertices are green, blue and red, respectively. 



neighbor pair of spins is ±3/2, and set h c — 10, which is 
larger than the largest possible dipolar field g acting on 
any island [13] . As long as h c > g, results with different 
h c are simply shifted relative to one another. 

We introduce the notation that © refers to a type i 
vertex, and rij the population fraction of type i vertices. 
In the bulk of the array, topological and energetic con- 
straints allow the following 2-vertex processes: 



@© -> ®®, ®@, (D®, ©@ 

as shown in Fig. [Tg. Because of the energetic constraints, 
type 4 vertices cannot be created from an initial polarised 
state, type 1 vertices cannot be destroyed, and the pro- 
cess ©© — !• ®® cannot occur. The impossibility of de- 
stroying type 1 vertices and the impossibility of the pro- 
cess ®® — > @® give rise to an important phenomenon, 
trapping, in which a type 3 vertex or a region of type 2 
vertices is frozen if surrounded by type 1 vertices. The 
existence of trapping contributes to a slowing of dynam- 
ics at long times. Type 3 vertices can only be nucleated 
and expelled at the array edges, via the processes ®^® 
and ©•(->■©, also shown in Fig. [Tg. 

Starting from an initial saturated configuration, the 
dynamics begin with the nucleation of type 3 vertices 
at the array edges perpendicular to the original polar- 
ization, via the process © — > ®. Then type 3 ver- 
tices can propagate via two processes, @© — > ®@ and 



), the former requiring a smaller applied field 
than the latter. 

In the absence of long-range interactions, there would 
be four field regimes: a very low field regime where the 
field is too small to effect a response from the spins, a very 
high field regime in which the magnetization tracks the 
applied field, and two nontrivial intermediate regimes. 
The difference between the two intermediate regimes is 
that in the first (the 'low field' regime), type 3 vertices 
can only propagate via the process @© — > ®® whereas 
in the second (the 'high field' regime), both type 3 propa- 
gation processes can occur. Long range interactions mod- 
ify this picture somewhat because the dipolar fields act- 
ing on spins in different locations are not equal, causing 
a spread in the applied field strength required for ver- 
tex processes. This leads to a more gradual crossover 
between the low and high field regimes. 

A full description of vertex dynamics should take into 
account correlations and long range interactions, but if 
vertices are sufficiently well-mixed then these should be- 
come negligible on average and a mean field description 
based on short range interactions should hold. We pro- 
pose a set of ordinary differential equations to describe 
an open edge array under these assumptions. Compar- 
ison with numerical simulations that take into account 
the full long range interactions (see below) shows that 
the mean field approach describes the crossover and high 
field regimes relatively well. 

The equations for the rii{9) are: 



hi = 4^13(1 - ^V/2)n 2F n 3F 
+ 2{2vii + v l2 ){l-jr/2)n\ F 
h 2F = -4z^i3(l - ,jV j2)n 2F n 3F 



vl V jV n 3F , 



2(v 12 + 2v 22 )(\ - J / /2)n 2 3F - v^-yYn 



L2F 



(la) 



(lb) 



+ f 32 ^Vn 3F — ^t(1 — ^V)n 2F n\, 
h 3F = -4(i/n + 1/12 + 1/22XI - ^/ 2 ) n 3F + v e 2Z ,jYn 2F 
~ O31 + vl 2 )jYn 3F - v T {\ ~ ,yV)n 3F nl, 

(1c) 

h iT = M 1 - ^)n lF n\ (t = 2, 3), (Id) 

where hi = dn^/dfl, jY is the fraction of vertices at the 
array edge, and subscripts 'T' and 'F' denote vertices 
that are trapped and vertices that are free to be involved 
in dynamics, respectively. The reaction rates v depend 
on the total field acting on the flipping spin. The dipolar 
part is calculated in the approximation that only spins 
belonging to the spin's two vertices contribute to it. 

We integrate the system (JTJ) with initial conditions 
n 2 p(0) — 1, n>i(0) = for all other i, for h = 10. In Fig. 
[2] we see that size effects are clearly significant for pop- 
ulation levels: as the number of vertices N increases the 
final type 1 and 3 populations decrease logarithmically 
slowly (see Inset), the number of rotations required to 
reach a steady state increases, and the transient peak in 
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FIG. 2: (a) Type 1 and (b) type 3 population fractions for 
dynamics given by Eqs. {TJ, f° r a field h = 10 and for arrays 
of increasing size. Inset: Final type 1 population fraction as 
a function of array size N. 



713 disappears. These effects can all be attributed to the 
decrease in JV with N, because for large arrays with rela- 
tively few edge vertices, the rate of type 3 vertex creation 
will be smaller, resulting in reduced rates for processes 
that create type 1 vertices. 

Considerations for closed edge arrays can be made in 
direct analogy to those presented above for open edge 
arrays. The resulting equations are more complicated 
due to the additional populations introduced by the edge 
geometry, as noted earlier, but the resulting behaviors 
for m and n% turn out to be very similar to those shown 
in Fig. [2] This similarity signals a failure of our mean 
field approximation. 

Correlations and long range dipolar effects. — The pri- 
mary shortcoming of our mean field population dynamics 
picture is the neglect of spin correlations, but the neglect 
of long range dipolar interactions is also important. We 
demonstrate this using numerical simulations which do 
not suffer from these two shortcomings. At each step 
of a simulation, all spins are accessed in a random order, 
and flipped according to the h c rule described above. The 
field is held constant until no further spin flips can occur 
and then rotated by d9 = 0.01. The initial configuration 
is saturated in the x direction. Example animated time 
evolutions for open and closed arrays can be found in the 
Supplementary Materials 14 j. 
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FIG. 3: (a) Type 1 and (b) type 3 mean vertex populations 
from numerical simulations, for h — 10.75 (solid black lines), 
10 (blue dashed lines), and 9.25 (red dotted lines), for an open 
edge array of 400 islands. Averages are over 100 runs of the 
simulation. The inset of (a) shows the dependence of the final 
type 1 population on h, with the two field regimes of interest: 
low fields (9 < h < 9.5) and high fields (10 < h < 11). 



In Fig. [3]we plot (ni(9)} and (n 3 (9)) for an open edge 
array. The inset shows the dependence of the final type 
1 population fraction on the field strength h. We find 
four field regimes, as expected by our general considera- 
tions of vertex dynamics. The lowest and highest regimes 
are trivial: for h < 9, no spin flips can occur, while 
for ft, > 11 spins simply follow the field. In the interval 
9 < ft. < 11, the low field regime (red dotted lines) and 
high field regime (solid black lines) are shown, together 
with a crossover between the two (dashed blue lines) . It 
is worth stressing we are able to attain (ni) values of up 
to 90% if the field is slightly larger than the lowest field 
allowing dynamics. 

We performed similar simulations for closed edge ar- 
rays. As an example of the results, Fig. 0] shows the 
evolution of mean n\ and 77,3 values for an applied field 
ft = 11.25, with an inset showing the dependence of the 
final (ni) on ft. These results are quite different to those 
obtained for open edge array for several reasons. 

The same 2- vertex processes occur in the bulk of open 
and closed arrays and likewise type 3 nucleation and ex- 
pulsion occurs at the edges of closed arrays. However, in 
closed edge arrays, type 3 nucleation and expulsion can 
only occur when the neighboring 3-island edge vertex is 
of type 2e, because the fields required when it is type le 
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FIG. 4: Type 1 (dotted line) and type 3 (solid line) vertex 
population fractions from simulations for a closed edge array 
of 220 islands with an applied field h = 11.25. Averages are 
over 100 runs of the simulation. Inset: Dependence of the 
mean final type 1 population fraction on field strength h. 



are prohibitively large. Thus, creation of type le vertices 
reduces the number of nucleation and expulsion sites and 
slows the dynamics. 

Moreover the fields required to nucleate type 3 vertices 
in the closed edge array are larger than the field required 
for the (D® — > ®@ process, so this process is always able 
to occur. As a result, there is only one interesting field 
regime, for 10.5 < h < 11.5. The field required for initial 
type 3 nucleation is sufficiently large that a relatively 
high proportion of nucleated type 3 vertices are driven 
across the array in the process @® — > ®@ and expelled 
at the opposite edge. This leads to oscillations in (7J3). 

Conclusions. — We have demonstrated in this paper 
that vertex dynamics can be formalized via population 
dynamics equations. We have done this within the lim- 
itations of a short range interaction mean field approxi- 
mation, and provided insights into array size dependence. 
Most significantly, we find that the final number of type 1 
vertices decreases very slowly with increasing array size. 

It may be possible to improve the mean field treat- 
ment by introducing spatially dependent densities for the 
populations and estimating effects of spin correlations 
on transition rates, but that is beyond the scope of the 
present work. However, we note that a mean field ap- 
proximation can provide an accurate description in some 
circumstances. The reason is the following. While evo- 
lution begins with nucleation of type 3 vertices at the 
array edge, numerically we find that different edges pro- 
duce very different results. Instead if the initial config- 
uration is random rather than saturated, many defects 
are included and edge effects become negligible. Thus a 
random initial configuration is described well in a mean 



field approach. 

Based on results from our models, we can also make a 
general observation of importance for the question of de- 
magnetization protocol. We find that a simple constant 
h protocol may be as good as more elaborate protocols, 
if h is correctly tuned. We have also studied the system 
using some protocols that change both 8 and h. The best 
demagnetization appears when h is varied such that the 
system spends a long time with h near the values that 
give the maximum mean n\ in the inset of Fig. [3^. 

Lastly, we also note that real islands may be modeled 
more accurately with Stoner-Wohlfarth switching. We 
have examined this and find that we obtain qualitatively 
similar results for critical field and Stoner-Wohlfarth 
switching, but critical field switching is more amenable 
to a simple theoretical description. 

Z.B. and R.L.S. thank the ARC and the WUN for sup- 
port. Z.B. also acknowledges ARNAM and the Hackett 
Foundation. 



[1 

[2 

[3; 

[4 

[5 

[6 
[7; 

[8 
[9 
[10 

[11 

[12 

[13 
[14 



M. J. Harris, S. T. Bramwell, D. F. McMorrow, T. Zeiske, 
and K. W. Godfrey, Phys. Rev. Lett. 79, 2554 (1997). 
S. T. Bramwell and M. J. P. Gingras, Science 294, 1495 
(2001). 

A. P. Ramirez, A. Hayashi, R. J. Cava, R. Siddharthan, 

and B. S. Shastry, Nature 399, 333 (1999). 

J. Snyder, B. G. Ueland, J. S. Slusky, H. Karunadasa, R. 

J. Cava, and P. Schiffer, Phys. Rev. B 69, 064414 (2004). 

S. V. Isakov, R. Moessner, and S. L. Sondhi, Phys. Rev. 

Lett. 95, 217201 (2005). 

R. F. Wang et al, Nature 439, 303 (2006). 

Y. Qi, T. Brintlinger, and J. Cumings, Phys. Rev. B 77, 

094418 (2008). 

J. Li, X. Ke, S. Zhang, D. Garand, C. Nisoli, P. Lammert, 
V. H. Crespi, and P. Schiffer, Phys. Rev. B 81, 092406 
(2010). 

X. Ke, J. Li, C. Nisoli, P. E. Lammert, W. McConville, 
R. F. Wang, V. H. Crespi, and P. Schiffer, Phys. Rev. 
Lett. 101, 037205 (2008). 

C. Nisoli, R. Wang, J. Li, W. F. McConville, P. E. Lam- 
mert, P. Schiffer, and V. H. Crespi, Phys. Rev. Lett. 98, 
217203 (2007). 

C. Niso li, J. Li, X. Ke, D . Garand, P. Schiffer, and V. H. 
Crespi, arXiv: 1005.3463 (unpublished). 
G. T. Barkema and M. E. J. Newman, Phys. Rev. E 57, 
1155 (1998). 

In the array studied by Ke et al. [jj, h c « 770 Oe and 
h dip » 10 Oe. 

Sec EPAPS Document X. For more information on 
EPAPS, see |http: / /www. aip.o rg/pubservs/epaps litml] 



